Figure 5 and 20: 1D plume evolution

Figure 5 and 20: 1D plume evolution#

For higher resolution figure, change npos in load-calc-hires-zoom-av.ipynb

%run import-modules-grid

import holoviews as hv
from holoviews import opts
hv.extension('bokeh')

from bokeh.plotting import figure, show
from bokeh.io import output_notebook

!jupyter --version

import os.path
import copy

output_notebook()
Platform:  Darwin Kernel Version 24.1.0: Thu Oct 10 21:00:32 PDT 2024; root:xnu-11215.41.3~2/RELEASE_ARM64_T6030
python version:  3.11.10
matplotlib version:  3.9.2
hvplot version:  0.11.0
numpy version:  2.1.2
pandas version:  2.2.3
pickle version:  4.0
yaml version:  6.0.2
pint version:  0.24.3
pyko version:  v0.8.3-dev-2024-05-12
print eos_table version:  v1.1.5b

Number of CPUs in the system: 12
Selected Jupyter core packages...
IPython          : 8.28.0
ipykernel        : 6.29.5
ipywidgets       : 8.1.5
jupyter_client   : 8.6.3
jupyter_core     : 5.7.2
jupyter_server   : 2.14.2
jupyterlab       : 4.2.5
nbclient         : 0.10.0
nbconvert        : 7.16.4
nbformat         : 5.10.4
notebook         : not installed
qtconsole        : not installed
traitlets        : 5.14.3
Loading BokehJS ...

Read in example vapor plume calculation#

rplumeinitarr = np.asarray([25.e3]) # m
rnebinitarr = np.asarray([250.e6]) #m
tsavearr = np.asarray([1.]) #s
pinitarr = np.flip(np.asarray([80.e9]))
velinitarr = np.asarray([0.]) # m/s

tstall = np.zeros((len(rplumeinitarr),len(pinitarr),len(velinitarr)))
if 0:
    fin='./data/temp/vp-h2o-tempgridlong.yml'
    fout='./data/temp/vp-h2o-tempgridlong-'

    ftype='YAML'
    verbose=True
    userdtstart=0.0001
    usertstepscale=1
    binoutput=True
    debug=False
    initarr=False
    #run = RunClass(fin=fin,fout=fout,ftype=ftype)    # initialize run parameters class
    #
    # read in the run template
    #readinput_yaml(run,verbose=0)
    #print(run.ieos[0].path)
    #hugarr = np.loadtxt(run.ieos[0].path+'NEW-SESAME-HUG-E.TXT',skiprows=3,delimiter=',')
    #print('Hugoniot file size = ',hugarr.shape)
    # vary starting parameters
    for ipp in range(len(pinitarr)):
        for irp in range(len(rplumeinitarr)):
            for ivel in range(len(velinitarr)):            
                fileid = 'P'+str(np.round(pinitarr[ipp]/1.e9))+'-R'+str(np.round(rplumeinitarr[irp]/1.e3))+'-V'+str(np.round(velinitarr[ivel]/1.e3))
                outputfilename = fout+fileid+'.dat'
                print("#outputfilename='"+outputfilename+"'")
                if os.path.isfile(outputfilename):
                    %run load-calc-hires-zoom-av.ipynb
                    print('loading '+outputfilename)
                else:
                    print('file not found: '+outputfilename)


    pkopostrack_orig = copy.deepcopy(pkopostrack)    
    # av version is to mitigate the artificial viscosity error at the material interface in the plot
    # note: material interface temperature error suppressed for plotting for +- 5 cells
    pkopostrack = copy.deepcopy(pkopostrack_av)

    outputfilename='./data/fig5-data.pkl'
    with open(outputfilename,"wb") as f:
        pickle.dump([pko,pkopostrack,rstall],f)
if 1:
    outputfilename='./data/fig5-data.pkl'
    with open(outputfilename,"rb") as f:
        [pko,pkopostrack,rstall] = pickle.load(f)
jmaxmat1=max(pko['j'][pko['time']==0][pko['mat']==1])
## this tracks the plume front versus time by following the index j at the edge

timeall = np.unique(pko['time'])
imat1 = np.where((pko['mat'] == 1))
jmaxmat1=max(pko['j'][pko['time']==0][pko['mat']==1])
jmat1all =np.where((pko['j']==jmaxmat1))[0]
posmat1all = (pko['pos'][jmat1all])
#print(jmat1all)
#plt.plot(timeall,posmat1all)
#print(timeall,posmat1all)
#edgeline = hv.Curve((timeall,posmat1all)).opts(color='white')
edgelinewtdashthin = hv.Curve((timeall,posmat1all)).opts(color='white',line_dash='dashed',line_width=.5)
edgelinewtdashthick = hv.Curve((timeall,posmat1all)).opts(color='white',line_dash='dashed',line_width=4)
#edgelinebl = hv.Curve((timeall,posmat1all)).opts(color='black')
edgelinebldashthick = hv.Curve((timeall,posmat1all)).opts(color='darkgrey',line_dash='dashed',line_width=4)
edgelinebldashthin = hv.Curve((timeall,posmat1all)).opts(color='darkgrey',line_dash='dashed',line_width=.5)
edgelinebldashthick
#edgelinebldash = hv.Curve((timeall,posmat1all)).opts(color='black',line_dash='dashed',line_width=.5)
print(max(posmat1all))
itmp = np.where(posmat1all == max(posmat1all))[0]
print(itmp)
print(timeall[itmp])
stall_t = timeall[itmp]
stall_p = max(posmat1all)
stall_pt = hv.Scatter((stall_t,stall_p)).opts(size=8, fill_color="grey", line_color="grey")
86.3635503798735
[428]
[11.87500197]
posline = 40*np.power(timeall,1)
possedov = 40*np.power(timeall,.4)
sedovline = hv.Curve((timeall,possedov)).opts(color='darkgrey',line_width=4)
sedovline2 = hv.Curve((timeall+.75,possedov)).opts(color='black',line_width=4)
startline = hv.Curve((timeall,posline)).opts(color='darkgrey',line_width=4,line_dash='dotted')
#edgelinebldashthick*sedovline*startline.opts(ylim=(0,100))
sedovplot=edgelinebldashthick*sedovline*sedovline2*startline*stall_pt.opts(ylim=(0,100),width=500,height=400,xlabel='Time (hr)',ylabel='Radial Position (Mm)',fontsize={'title': 16, 'labels': 14, 'xticks': 12, 'yticks': 12, 'cticks':12})
#sedovplot=edgelinebldashthick*sedovline*sedovline2*startline*stall_pt.opts(ylim=(0,100),width=1200,height=800,xlabel='Time (hr)',ylabel='Radial Position (Mm)',fontsize={'title': 26, 'labels': 24, 'xticks': 22, 'yticks': 22, 'cticks':22})
sedovplot
# Figure 20 was saved as a png from the widget below
hv.save(sedovplot,'./plots/Fig20-sedovplot.png')
# for each initial position, track
#  0 - time of first nebula shock
#  1 - time of first plume front
#  2 - peak velocity nebula shock and time
#  3 - peak velocity of plume shock and time
#  4 - peak temperature of nebula shock and time (not the same as velocity)
#  5 - peak temperature of plume shock and time (maybe same as velocity)
#  6 - peak pressure of nebula shock and time (not the same as velocity)
#  7 - peak pressure of plume shock and time (maybe same as velocity)
#  8 - time with first shock

npos = 1024
waveinfoarr = np.zeros([9,npos])
waveinfoarr[6,:]=np.log10(1.4064186900465516E-07)
tind=[]
tind2=[]


tracktimeall = np.unique(pkopostrack['time'])
ntime = len(tracktimeall)
itimestall = np.where(tracktimeall > tstall/3600.)[0][0]
trackposall = np.unique(pkopostrack['pos'])
iposstall = np.where(trackposall < rstall/1.e6)[0]
for ipos in np.arange(2,len(trackposall)):
    tind=[]
    tind2=[]    
    # first motion?
    pkoposone = pkopostrack[pkopostrack['pos']==trackposall[ipos]]
    pkoposonemat1 = pkoposone[pkoposone['mat']<1.5]
    pkoposonemat2 = pkoposone[pkoposone['mat']>1.5]
    ind1= np.where(np.asarray(pkoposone['up']>10))[0] # moving positions
    if len(ind1)>0:
        waveinfoarr[0,ipos]=tracktimeall[ind1[0]]
        #waveinfoarr[8,ipos]=
        #print(tracktimeall[ind1[0]])

    ind2 = np.where(np.asarray(pkoposone['mat']==1))[0] # mat1 positions
    if len(ind2)>0:
        waveinfoarr[1,ipos]=tracktimeall[ind2[0]]
        #print(ind1[0],ind2[0])
        if len(ind1)>0:
            tind = np.arange(ind1[0],ind2[0]-1)
        tind2 = np.arange(ind2[0],itimestall)
    else:
        if len(ind1)>0:
            tind = np.arange(ntime)
        else:
            tind=[]#np.arange(npos)

    if len(tind)>0:
        #waveinfoarr[2,ipos] = np.max(np.asarray(pkoposone['up'])[tind])
        #waveinfoarr[4,ipos] = np.max(np.asarray(pkoposone['temp'])[tind])
        waveinfoarr[2,ipos] = np.max(pkoposonemat2['up'])
        waveinfoarr[4,ipos] = np.max(pkoposonemat2['temp'])
        waveinfoarr[6,ipos] = np.max(pkoposonemat2['rho'])
    if len(tind2)>0:
        #waveinfoarr[3,ipos] = np.max(np.asarray(pkoposone['up'])[tind2])
        #waveinfoarr[5,ipos] = np.max(np.asarray(pkoposone['temp'])[tind2])
        waveinfoarr[3,ipos] = np.max(pkoposonemat1['up'])
        waveinfoarr[5,ipos] = np.max(pkoposonemat1['temp'])
        waveinfoarr[7,ipos] = np.max(pkoposonemat1['rho'])
s1len = len(np.where(timeall<17)[0])
s1posall = np.zeros([s1len])
for it in np.arange(s1len):
    temp = pkopostrack[pkopostrack['time']==timeall[it]]
    temp2 = temp[temp['temp']>250]
    #print(it,min(temp2['pos']))
    s1posall[it]=(min(temp2['pos']))

#    display(temp2)
#    print(temp[pkopostrack['temp']>400][0])

#plt.plot(timeall[0:s1len],s1posall)
s2frontblthick = hv.Curve((timeall[15:s1len],s1posall[15::])).opts(color='black',line_dash='dotted',line_width=4)
s2frontwtthick = hv.Curve((timeall[15:s1len],s1posall[15::])).opts(color='black',line_dash='dotted',line_width=4)
s2frontblthick
#s1front = hv.Curve((timeall,waveinfoarr[2,:])
print(len(s1posall))
print(len(timeall[0:s1len]))
505
505
s1len = len(np.where(timeall<6.3)[0])
s1posall = np.zeros([s1len])
for it in range(s1len):
    temp = pkopostrack[pkopostrack['time']==timeall[it]]
    temp2 = temp[temp['up']>0.8]['pos']
    if len(temp2)>1:
        s1posall[it]=max(temp2)
    

#    display(temp2)
#    print(temp[pkopostrack['temp']>400][0])

#plt.plot(timeall[0:s1len],s1posall)
s1frontblthick = hv.Curve((timeall[0:s1len],s1posall)).opts(color='black',line_dash='solid',line_width=4)
s1frontwtthick = hv.Curve((timeall[0:s1len],s1posall)).opts(color='white',line_dash='solid',line_width=4)
s1frontblthick
#s1front = hv.Curve((timeall,waveinfoarr[2,:])
## this tracks the plume front versus time by following the index j at the edge

jradstall2=min(pko['j'][pko['time']==0][pko['pos']>rstall/1.e6*1.25])
print(jradstall2)
streamline0 = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='darkgrey',line_dash='solid',line_width=1.5)
streamline0wt = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='white',line_dash='solid',line_width=1.)


jradstall2=min(pko['j'][pko['time']==0][pko['pos']>rstall/1.e6])
print(jradstall2)
streamline1 = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='darkgrey',line_dash='solid',line_width=1.5)
streamline1wt = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='white',line_dash='solid',line_width=1.)

jradstall2=min(pko['j'][pko['time']==0][pko['pos']>rstall/1.e6*.75])
print(jradstall2)
streamline2 = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='darkgrey',line_dash='solid',line_width=1.5)
streamline2wt = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='white',line_dash='solid',line_width=1.)

jradstall2=min(pko['j'][pko['time']==0][pko['pos']>rstall/1.e6*.5])
print(jradstall2)
streamline3 = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='darkgrey',line_dash='solid',line_width=1.5)
streamline3wt = hv.Curve((pko['time'][pko['j']==jradstall2],pko['pos'][pko['j']==jradstall2])).opts(color='white',line_dash='solid',line_width=1.)

edgelinebldashthick*streamline1*streamline2*streamline3*streamline0
543
475
407
337
# bigger fonts

# extract data from pko dataframe for an x-t diagram


fs1=18 # title
fs2=18 # labels
fs3=18 #xticks,yticks
fs4=16 #cticks

xrr=(0,14)
yrr=(0,100)

wpix = 800
hpix = 500

heatmap1 = pkopostrack.hvplot.heatmap(y='pos', x='time', C='up',
    title="A. Velocity Evolution", 
    clabel='Radial velocity (km/s)',clim=(-4,4),
    xlabel='Time (hr)',xlim=xrr,ylim=yrr,
    ylabel='Radial Position (Mm)',
    width=wpix, height=hpix).opts(cmap='bwy',fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})
    
heatmap1b = pkopostrack.hvplot.heatmap(y='pos', x='time', C='up',
    title="A. Vapor Plume Features", 
    clabel='Radial velocity (km/s)',clim=(-4,4),
    xlabel='Time (hr)',xlim=xrr,ylim=yrr,
    ylabel='Radial Position (Mm)',
    width=wpix, height=hpix).opts(cmap='bwy', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})

heatmap2 = pkopostrack.hvplot.heatmap(y='pos', x='time', C='pres',
    title="B. Pressure Evolution", 
    clabel='log₁₀P (Pa)',
    xlabel='Time (hr)',
    ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
    clim=(-5,1),
    width=wpix, height=hpix).opts(cmap='bgy', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})

heatmap2b = pkopostrack.hvplot.heatmap(y='pos', x='time', C='pres',
    title="B. Distinct Zones Followed by Mixing", 
    clabel='log₁₀P (Pa)',
    xlabel='Time (hr)',
    ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
    clim=(-5,1),
    width=wpix, height=hpix).opts(cmap='bgy', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})

heatmap3 = pkopostrack.hvplot.heatmap(y='pos', x='time', C='rho',
    title="C. Density Evolution", 
    clabel='log₁₀ρ (g/cm³)',
    xlabel='Time (hr)',
    ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
    width=wpix, height=hpix,clim=(-9,-5)).opts(cmap='bgy', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})

heatmap3b = pkopostrack.hvplot.heatmap(y='pos', x='time', C='rho',
    title="C. Size Sorting", 
    clabel='log₁₀ρ (g/cm³)',
    xlabel='Time (hr)',
    ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
    width=wpix, height=hpix,clim=(-9,-5)).opts(cmap='bgy', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4},xlim=(0,13))

heatmap4 = pkopostrack.hvplot.heatmap(y='pos', x='time', C='temp',
    title="D. Temperature Evolution", 
    clabel='T (K)',
    xlabel='Time (hr)',
    ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
    width=wpix, height=hpix,clim=(200,2000)).opts(cmap='plasma', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})

heatmap4b = pkopostrack.hvplot.heatmap(y='pos', x='time', C='temp',
    title="D. Thermal and Redox Zones", 
    clabel='T (K)',
    xlabel='Time (hr)',
    ylabel='Radial Position (Mm)',xlim=xrr,ylim=yrr,
    width=wpix, height=hpix,clim=(200,2000)).opts(cmap='plasma', fontsize={'title': fs1, 'labels': fs2, 'xlabel': fs3, 'ylabel': fs3, 'cticks':fs4,'xticks':fs4,'yticks':fs4})



linegroup=edgelinebldashthick*stall_pt*s2frontblthick*s1frontblthick
linegroupwt=edgelinewtdashthick*stall_pt*s2frontwtthick*s1frontblthick
streamlines=streamline1*streamline2*streamline3*streamline0
streamlineswt=streamline1wt*streamline2wt*streamline3wt*streamline0wt
#heatmap*edgeline
#group=(heatmap*edgelinedash+heatmap1*edgelinebldash+heatmap2*edgelinedash+heatmap3*edgelinedash).cols(2)
grouplabels2=(heatmap1b*linegroup*streamlines+heatmap2b*linegroupwt*streamlineswt+heatmap3b*linegroupwt*streamlineswt + heatmap4b*linegroupwt*streamlineswt).cols(2)

grouplabels2
hv.save(grouplabels2, './plots/Fig5-1dsim.png')

Interactive version of Figure 5#

The code below generates the html interactive with the results from this 1D simulation: link to online interactive version.

#xrr=[0,150000]
colstr = [ 'blue','orange','grey']
xrr=[0,250]
densityplot = pko.hvplot.scatter(x='pos',y='rho',groupby='time',by='material',hover_cols=['all'],
                       #title='pyKO vs Time (seconds)',
                       logy=True,group_label='test',
                       xlabel='Radial Position (Mm)',ylabel='Density (g/cm³)',xlim=xrr,ylim=[1.e-12,10],color=colstr,fontsize={'title': 16, 'labels': 14, 'xticks': 12, 'yticks': 12, 'legend':12,'legend_title':14})

tempplot = pko.hvplot.scatter(x='pos',y='temp',groupby='time',by='material',hover_cols=['all'],
                      # title='pyKO vs Time (seconds)',
                              xlabel='Radial Position (Mm)',ylabel='Temperature (K)',xlim=xrr,ylim=[0,4000],color=colstr,fontsize={'title': 16, 'labels': 14, 'xticks': 12, 'yticks': 12, 'legend':12,'legend_title':14})

upplot = pko.hvplot.scatter(x='pos',y='up',groupby='time',by='material',hover_cols=['all'],
                       #title='pyKO vs Time (seconds)',
                    title='Water plume R=25 km, P=80 GPa into 0.1 Pa',

                       xlabel='Radial Position (Mm)',ylabel='Radial Velocity (m/s)',xlim=xrr,ylim=[-3000,10000],color=colstr,fontsize={'title': 16, 'labels': 14, 'xticks': 12, 'yticks': 12, 'legend':12,'legend_title':14})

pressureplot = pko.hvplot.scatter(x='pos',y='pres',groupby='time',by='material',
                       title='Time in hours',
                       xlabel='Radial Position (Mm)',ylabel='Pressure (Pa)',logy=True,ylim=[1.e-6,1.e11],xlim=xrr,color=colstr,fontsize={'title': 16, 'labels': 14, 'xticks': 12, 'yticks': 12, 'legend':12,'legend_title':14})

hv.output(widget_location='top')
groupplot = (upplot+ pressureplot + densityplot+tempplot ).cols(2)
groupplot

Material 1 = water plume
Material 2 = nebular gas within the stall radius
Material 3 = nebular gas beyond the stall radius

# save the interactive html file
# this save takes a while so comment out unless needed
# hv.save(groupplot, './plots/Fig5-interactive.html')